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We describe an algorithm and a C++ implementation that we have written and made available for 
calculating the fully nonlinear evolution of 5D braneworld models with scalar fields. Bulk fields allow 
for the stabilization of the extra dimension. However, they complicate the dynamics of the system, 
so that analytic calculations (performed within an effective 4D theory) are usually only reliable for 
static bulk configurations or when the evolution of the extra dimension is negligible. In the general 
case, the nonlinear 5D dynamics can be studied numerically, and the algorithm and code we describe 
are the first ones of that type designed for this task. The program and its full documentation are 
available on the Web at http://www.cita.utoronto.ca/~jmartin/BRANECODE/ 1 . In this paper we 
provide a brief overview of what the program does and how to use it. 



I. INTRODUCTION 

Many extensions of the Standard Model have in com- 
mon the presence of extra dimensions. This has to be 
contrasted with the fact that our world looks four di- 
mensional, so one has to explain why the presence of the 
extra space has not yet been detected. The traditional 
answer has been that the extra space is compact and very 
small, so that the fields associated with its excitations are 
too heavy to be observable in accelerators or cosmology. 
More recently, it has been realized that ordinary matter 
and gauge interactions may be confined on lower dimen- 
sional submanifolds, known as branes. In this case, they 
could be four dimensional objects, even if the geome- 
try of the theory is higher dimensional. The situation 
is different for gravity, which propagates in the whole 
bulk space. Several questions naturally arise, such as 
why a compact space would remain small while the three 
non-compact dimensions are undergoing cosmological ex- 
pansion, or why the expansion of the universe we see is 
described by 3 + 1 dimensional general relativity so well. 
The presence of extra dimensions may cause deviations 
from the standard FRW cosmology that is supported by 
observations. 

In most cases, these two questions turn out to be inti- 
mately related. Only if the extra space is static can the 
evolution of the non-compact coordinates behave as in 
the standard four dimensional case. Hence, the dynam- 
ics of the hidden dimensions becomes a crucial ingredient 
in understanding the evolution of the ones we observe. In 
some particular cases, static bulk configurations can be 
achieved under the combined action of the bulk/brane 
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gravity. In most realistic examples that could account 
for our observed four dimensional cosmology the stabil- 
ity is due to the presence of additional fields that acquire 
nontrivial configurations in the bulk. While the stabiliza- 
tion has to be effective at relatively "late" times, the first 
stages of our universe (before primordial nucleosynthesis, 
for instance) are much less constrained. The evolution of 
the bulk may have been significant at this phase, and this 
offers many new possibilities for phenomenology. This is 
particularly true with the addition of the fields responsi- 
ble for the "late" time stabilization, since they constitute 
new dynamical degrees of freedom for the system. 

While the above considerations are valid for all mod- 
els with extra dimensions, significant computations have 
been performed in the framework of brane models. These 
models can be thought of as simplified, phenomenologi- 
cal (bottom-up) versions of branes in string theory. The 
string dynamics is ignored and the primary focus is on 
the classical dynamics from the point of view of Gen- 
eral Relativity. Branes act as a source for the Einstein 
equations of the system, with their tension and possibly 
with the energy density of fields confined on them. Ad- 
ditional sources are a bulk cosmological constant, or the 
energy density of possible bulk fields. This set-up is suf- 
ficiently rich to describe very interesting situations. For 
example, inflation in braneworlds can acquire a nice geo- 
metrical interpretation, with the inflaton associated with 
the distance between different branes, the Hubble param- 
eter scale associated with the induced curvature on the 
branes, and with reheating through radion oscillations. 

Despite this great simplification, the whole dynamics 
is still very complicated, in particular when a bulk scalar 
field is present. In this case, analytical computations are 
typically performed within an effective 4 dimensional the- 
ory, obtained after the extra dimension is integrated out, 
or perturbatively, using linearized analysis around simple 
backgrounds. While these studies are very useful when 
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the extra space is static or quasi-static, they are not suf- 
ficient to describe the system when the evolution of the 
bulk is important. For example, systems that are stable 
at low energy (low curvature of the branes) can become 
unstable when the energy/curvature is increased. The 
stability/instability can be studied analytically. How- 
ever, one cannot determine analytically where the system 
will evolve towards when the initial configuration turns 
out to be unstable. 

For example, we numerically examined the dynamics 
of brane collisions and found that, as the branes approach 
each other, the spacetime of the bulk asymptotically ap- 
proaches the Kasner-type solution. 

Motivated by these limitations with the analytical 
treatment, we undertook a numerical study of these mod- 
els. We developed a numerical algorithm, implemented 
in a C++ code, specifically designed for codimension 
one braneworlds, with a scalar field included in order 
to provide stabilization of the bulk at low energies. With 
the assumption of homogeneity and isotropy along the 
brane spatial coordinates x (corresponding to the stan- 
dard assumption of homogeneity and isotropy of the non- 
compact coordinates) , the problem is reduced to an effec- 
tively 2 dimensional one. The independent variables are 
the bulk time t and the bulk dimension y. The program 
integrates numerically the full set of Einstein equations 
in the bulk, together with the Israel junction conditions 
at two orbifold branes. We discretize the two dimensional 
spacetime (time and the bulk coordinate) and solve the 
bulk differential equations by finite-differencing them us- 
ing the so-called leapfrog scheme. This algorithm is suffi- 
cient for our problem, and it provides a good compromise 
between accuracy and computational time (see Section 
[H]for more details). The two branes act as (one dimen- 
sional) boundaries of this space, and the junction condi- 
tions provide the boundary conditions for the system at 
each time-step. The solution of a boundary value prob- 
lem is required to provide generic initial conditions that 
fulfill the constraint equations and the boundary condi- 
tions at the beginning of the evolution. In the static 
configurations we have mentioned above, this boundary 
value problem is significantly simpler than in the general 
case. The setup of initial conditions is explained in more 
detail in Section HVI 

The first results obtained with the code have 
been presented in [|. We are now making 
the code public on the World Wide Web un- 
der the name BRANECODE. Its website is at 
http : //www. cita.utoronto . ca/~kof man/BRANECDDE/. 
The website for the program has documentation, includ- 
ing derivations of all the equations used in the program. 
Here we present a short summary of what the program 
does and what it can be used for. For more details see 
the website. Section [H] of this paper gives an overview of 
what the program is and how it works. Sections II I II and 
IIVI describe the evolution equations and the setting of 
initial conditions respectively. Section Ivl describes some 
of the output generated by the program. The references 
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II. OVERVIEW AND USER ADJUSTABLE 
FILES 

In this Section we give an overview of the program and 
how to adjust it for a particular simulation. More details 
can be found in the documentation. 

To work with the BRANECODE the user must specify 
a model, consisting of bulk and brane potentials for the 
scalar field, plus initial conditions for the field and the 
geometry. This information is encoded in a model file, 
which is a header file read in by the program. The file 
should be called model_numeric .h or model_analytic .h 
depending on whether the initial conditions are speci- 
fied numerically or analytically. The model files con- 
tain the potentials and their first and second deriva- 
tives that are needed for the evolution of the bulk equa- 
tions Q and boundary conditions JSJ. For example, the 
BRANECODE distribution includes both a numeric and 
an analytic model file (with different initial conditions). 
The file model_analytic . h contains examples for branes 
in AdS and AdS-Schwarzschild geometries, whereas the 
file model_aumeric .h we designed for the class of models 
with bulk scalar fields determined by the bulk and brane 
potentials 

V{4>) = \ m 2 4> 2 + A , 

U^) = iMi(& -o-if + Xi . (1) 

A bulk cosmological constant and brane tensions is in- 
cluded as constant terms in these potentials. Most im- 
portantly, the program is designed to work for arbitrary 
potentials, different from JQ. Other potentials V(<p), 
Ui{4>) can be implemented by modifying the correspond- 
ing lines in the model file. Aside from this file, the only 
other file that the user needs to modify is parameters .h, 
which contains all the parameters for a given run of the 
program. These include the number of grid points, the 
running time, and a number of other general variables 
specific to each run. There is also a parameter in this file 
that tells the program which type of model file to look 
for. 

Given a specific model and set of parameters, the 
BRANECODE solves the system of equations of motions 
for the metric functions and the scalar field Q along with 
the boundary conditions provided by the presence of the 
branes 10. The required functions are contained in the 
file equations . cpp. 

The BRANECODE has built-in routines for out- 
putting and plotting the metric fields, the scalar field, 
and derived quantities. These outputs are stored in 
ASCII files that can be read in and plotted by any 
standard plotting software. There are also options (set 
in parameters . h) to have the program call GNUPLDT 
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to generate and display postscript plots of the data at 
runtime. (These options should only be chosen if the 
program is running on a computer with both GNUPLDT 
and GHOSTVIEW.) An example of this graphical output is 
shown in Fig. ^ 
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--e^V ,(3) 
-e 2S ^ , 



B(y,t) 



supplemented by the constraint equations 




FIG. 1: An example of the graphical output generated by 
the BRANECODE. The plot shows the evolution of a met- 
ric component B (describing the interbrane separation) as a 
function of space and time. Physically, this evolution shows 
a transition from an unstable static warped geometry solu- 
tion towards a stable static solution. During this non-linear 
reconfiguration the interbrane distance and the Hubble scale 
of the de Sitter geometry decreases. 

Once all parameters have been set and you have mod- 
ified or created a model file according to your wishes you 
simply compile and run the BRANECODE. The code is 
designed to be platform independent and should work 
with any CH — h compiler. The makefile that comes with 
the distribution has entries and flags for the GNU gcc 
compiler and the INTEL ice compiler. You can select 
one of them or edit the makefile to invoke your favorite 
compiler. 



III. ALGORITHM 

The evolution equations solved by the BRANECODE 
are the set of Einstein/scalar field equations on an effec- 
tively 2 dimensional spacetime obtained after imposing 
homogeneity and isotropy on the non-compact spatial co- 
ordinates x . In [l| we showed that under these conditions 
it is always possible to choose coordinates that bring the 
metric to the form 
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with the two branes fixed at y = 0, 1 , respectively. This 
choice is motivated by the fact that in this coordinate sys- 
tem the lattice size is time independent due to the fixed 
position of the branes and the equations simplify signif- 
icantly. In this gauge, we have the following dynamical 



-A'A + B'A + A'B-A' = ~ 
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2A' 2 - A'B' + A" -A 2 -AB = -~4> 2 - ~4>' 2 - - e 2B V . 
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Dots and primes denote derivatives with respect to the 
time t and the coordinate y along the bulk, respectively. 
In addition, the program imposes the following junction 
(Israel) conditions at the positions of the branes 



A' = --Ue B 
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B' = --Ue B 
6 



(5) 



These junction conditions are equivalent to extending the 
space beyond the two branes, and imposing Z2 symmetry 
across the branes. 

In the coordinate system the characteristic prop- 
agation speed of the dynamical equations is always 1. 
This is advantageous from a numerical viewpoint as the 
size of the time step can be optimized uniformly by set- 
ting it equal to spatial grid separation. We describe 
below our implementation of the leapfrog discretization 
scheme, which is stable, second order accurate, and non- 
dissipative. 

The program discretizes the 2 dimensional {t, y} space 
and computes the value of the three functions B, A, <f> at 
each grid point. For any fixed time, the grid is made of 
N + 1 points equally spaced along the bulk, with and 
N corresponding to the locations of the two branes. The 
value of TV is set in parameters . h. The same grid spacing 
is taken in the y and t directions. The initial conditions 

are in the form \A(y),B(y), <j)(y),A(y), B(y), <p(y)j , and 

they can be chosen arbitrarily subject to the fact that 
they satisfy (@J and (JSJ) at the initial time. Of particular 
interest are configurations with an initially static bulk, 
since the program can be used to verify their stability 
and to study the dynamics when they are unstable as 
described in section llVl 

The discretized initial data (^Ai , Ai, Bi, Bi, (f>i , 

4>i*j , i = 0, . . . , N has to be converted to a form suit- 
able for the leapfrog algorithm. Instead of having the 
fields and their velocities at the time t = to, the algo- 
rithm needs the spatial profiles of the fields (Ai, Bi, <fii) 
at two subsequent moments in time t — to and t = to + e. 
In some special cases the profiles at the two initial times 
can be calculated analytically, but in general the program 
takes a Runge-Kutta step using the initial derivatives to 
calculate the field values at to + e. 
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The field values at all subsequent time steps are calcu- 
lated as follows. (Note that we will use / here to denote 
a generic field, i.e. for equations that apply to all three 
fields A, B, and <j>.) In the gauge chosen, the bulk equa- 
tions only contain derivatives in the form / — / , and 
f 9~ f 9' ■ Recall that our grid spacing e = 1/N is equal 
to our time step. Thus, to second order accuracy, we 
can write the derivatives at a given point in spacetime in 
terms of the values of its neighbors as 

/-/" = \ (/up + /d„ - fit - /rt) + O (e 2 ) , 



f <J - ./"// : : ^2 [^ U P ~~ ^ dn ) ^np ~ 9dn) + 

(/rt - fit) (9rt - git)} + O (e 2 ) , (6) 



where the indices label relative grid positions as defined 
in Fig. 12 In this way, the three differential equations (0} 
become three algebraic equations for the three unknown 
quantities {B up , ^4 up , (f> up ) . Notice that this can be done 
independently site by site in the bulk. 

Once the bulk field values at the new time have been 
determined, we can apply the junction conditions to 
advance the boundary field values. We discuss only the 
computation at the first brane (i = ) . The treatment 
for the second brane is analogous. Only first derivatives 
in y enter in (j3J). To preserve second order accuracy, we 
use an "asymmetric" discretization of the derivatives 



fo 



1 



(-3/ + 4/ 1 -/ 2 ) + O(e 2 ) 



(7) 



where subscripts indicate the grid position of each field 
value (see Fig. EJ. The junction conditions © thus be- 
come a set of three algebraic equations in terms of three 
unknowns. We can eliminate Aq and Bq in favor of (f>o , 
and write an equation in terms of the only unknown 
quantity <j)Q , 



4^-02-300 -« Bl(W Wo) = 



(8) 



For specific brane potentials, these equations can be 
solved analytically. However, as the code is designed 
for arbitrary bulk/brane potentials, eq. (jHJ is solved nu- 
merically using the iterative Newton's method. Once 0o 
is determined, the remaining unknowns Bq and Aq are 
trivially computed through the remaining junction con- 
ditions. 



IV. INITIAL CONDITIONS 

In general, the specification of initial conditions, i.e. 
the determination of the initial spatial profiles for the 
fields and their velocities A(y), B(y), <j)(y), A(y), B(y), 
and 4>(y) is a non-trivial task. These function cannot be 
chosen independently, but rather are subject to the con- 
straint equations Q and boundary conditions (j3J) ■ How- 
ever, thephysical instability of static de Sitter configura- 
tions provides the possibility to generate interesting 
braneworld dynamics with static solutions as initial con- 
ditions that we can recommend. All the static solutions 
for a given model can be exhaustively classified by the 
phasespace analysis of the dynamics of the gravity /scalar 
system performed in 0. 

One simple way to generate static initial conditions is 
to consider a configuration of the form B — B (y) , <j) = 
4>(y) , A = B + H t (de Sitter branes in a static bulk) . 
The first of the equations Q is then trivially satisfied. 
For a sufficiently simple model, the remaining constraint 
equation can be solved analytically. Otherwise we pro- 
vide a MATHEMATICA notebook and MAPLE worksheet to 
solve them numerically for a given bulk potential and 
generate appropriately formatted initial data. For sim- 
plicity, the code determines two of the three parameters 
of the brane potentials Q in such a way that the bound- 
ary conditions (JSJ) are fulfilled initially. The third param- 
eter, either the tension on the brane Xi or the minimum 
of the brane potential <7j, is set by the user in the file 
model_numeric .h. If the user instead wants to set up 
initial conditions for a given choice of brane potentials, 
he can make use of a shooting method operating in the 
phasespace of the static solutions. Namely, he can start 
with initial conditions chosen such that they fulfill the 
boundary conditions on one of the branes, compute the 
corresponding bulk configuration, and keep varying the 
initial conditions (e.g. with a numerical scan) until the so- 
lution also satisfies the junction conditions at the second 
boundary. For any given set of bulk and brane poten- 
tials, there can be none, one, or more than one solution 
to the boundary value problem. The latter case opens 
up the possibility for interesting dynamics of transitions 
as investigated in p|. 
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FIG. 2: Numerical evolution scheme. See text for details. 



OUTPUT 



The main outputs of the program are the values of 
the three functions B, A, and <fi at different bulk sites 
and time steps. Two parameters inside parameters. h 
control how many points (both in the y and t directions) 
are to be saved. From these quantities, one can construct 
some outputs of immediate physical relevance. One is 
the physical interbrane distance. The branes are at a 
fixed coordinate distance in our gauge and their physical 
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separation is encoded in the metric coefficient B . 



D{t) 



dye 



B(t,y) 



(9) 



Another interesting quantity is the Hubble parameter as 
computed by observers on each of the two branes 



1 da 

tii = - — 

a dr 



At 



(10) 



where r is the physical time measured by the observers, 
defined by Dt = e B dt . Quite interestingly, the gauge 
choice used in the algorithm (see the previous Sec- 
tion) does not exhaust the gauge freedom of the prob- 
lem. Residual gauge transformations have been described 
in , and one can show that they do not change the val- 
ues of D and Hi , which therefore have physical meaning. 
Besides these physical quantities the BRANECODE also 
computes the Ricci scalar and the square of the Weyl 
tensor in the bulk. Other quantities of interest can easily 
be obtained from the "raw" values of B,A,d>. 



VI. CONCLUSIONS 

The main motivation for our work was to extend the 
knowledge of braneworld dynamics beyond the few situa- 
tions where it was known analytically. Apart from these 
situations (characterized by a static or slowly evolving 
bulk), approximate methods based on effective 4 dimen- 
sional computations are unreliable. In this short note 
we have presented an algorithm, together with its CH — h 
implementation, designed for numerical computations in 
this framework. First results obtained from this code 
were presented in 0, Q . We could show that some bulk 
configurations which are stable at low energy (low value 
of the expansion rate H of the two branes) become un- 
stable as H increases, in agreement with the analytical 
calculations of |3|]. This can be interpreted as a part of 
a more general phenomenon of gravitational instability 



of compactification to four dimensional de Sitter geom- 
etry yl- The numerical integration allowed us to fol- 
low the evolution of the system starting from the un- 
stable configuration. For certain bulk/brane potentials 
the system may evolve towards another static, but stable 
configuration, characterized by a lower value of H. The 
transition is typically a process of quick bulk reconfigu- 
ration. In many cases, however, the second configuration 
does not exist or cannot be reached. The two branes 
then either move apart to infinite distance, or they col- 
lide. Brane collisions is a very interesting subject by 
itself. The study of the numerical evolution led us to the 
conclusion that the asymptotic geometry is given by the 
universal Kasner-type that is typical of homogeneous but 
anisotropic strong gravity regimes. We did not investi- 
gated the case of branes departing from each other when 
one of them in general approaches a naked singularity in 
the bulk configuration. We also did not studied the pos- 
sibility of the formation of an apparent horizon between 
the branes. 



Our algorithm is focused on the simplest possible set- 
up which allows for brane stabilization based on a gener- 
alized Goldberger-Wise mechanism. Only one bulk scalar 
field has been considered, although with arbitrary poten- 
tials in the bulk and at the two branes. The inclusion of 
more scalar fields would be straightforward. In particu- 
lar, one could consider other fields, which are confined to 
the branes, and which are coupled to the bulk fields e.g. 
through the brane potentials. (The interplay between 
bulk/brane fields may lead to novel interesting features 
not considered in Q). Another easy generalization would 
be the inclusion of perfect fluids on the branes (for ex- 
ample, describing standard matter and radiation). Less 
trivial but more interesting extensions could be the inclu- 
sion of other types of fields (form fields in the bulk, for 
instance), or evolution with more dimensions included. 
For example, relaxing the hypothesis of homogeneity and 
isotropy along the ordinary dimensions would allow the 
study of inhomogencous perturbations of the system. 
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